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Summary 

A new algorithm is introduced to compute finite- 
amplitude states using primitive variables for Rayleigh- 
Benard convection on relatively coarse meshes. The algo- 
rithm is based on a finite-difference matrix-splitting 
approach that separates all physical and dimensional 
effects into one-dimensional subsets. The nonlinear pat- 
tern selection process for steady convection in an air-filled 
square cavity with insulated side walls is investigated for 
Rayleigh numbers up to 20,000. The heated lower 
boundary is augmented with “noisy” boundary conditions 
to illustrate the transient and bimodal nature of the pattern 
selection process above the first critical Rayleigh number. 
Above a second critical Rayleigh number other instability 
modes may also be excited. The internalization of 
disturbances that evolve into coherent patterns is 
investigated and transient solutions from linear 
perturbation theory are compared with and contrasted to 
the full numerical simulations. 

The basins of attraction for transient solutions in a phase 
plane are shown to be bounded by singular states repre- 
senting unstable modes. A particular symmetric mode that 
may be stable to symmetric or random disturbances but 
unstable to antisymmetric disturbances is investigated. 


Introduction 

Rayleigh-Benard convection is considered by many an 
ideal model for hydrodynamic stability since convection 
and dissipation compete with one another in a relatively 
straightforward manner. Among the many proposed 
approximations to the Navier-Stokes equations, the sim- 
plified form first advocated by Boussinesq as a model of 
natural convection is now commonly used as the basis for 
most theoretical and numerical stability research. The 
classical solution method of linearization and expansion 
in sinusoidal functions as described in detail by 
Chandrasekhar (1961) is a firmly established tool for vir- 
tually all stability problems. The eigenvalue problem that 
determines the normal modes and critical parameters (in 
this case the Rayleigh number) is itself quite a difficult 
mathematical problem. In the continual search for simpler 
formulations, series expansions in terms of orthogonal 
functions other than sinusoids have been used with some 
success. In a recent study, Lee et al. (1989) used an 
expansion in Chebyshev polynomials to consider the 
linear stability of Rayleigh-Benard convection in a general 
two-dimensional cavity. The first instability mode is a 
single roll creating a vortex-type motion for cavities with 
aspect ratio less than two. The computed critical Rayleigh 
number for these finite-aspect-ratio cavities is well above 
the value 1708, which is appropriate for infinite-aspect- 


ratio systems. The normal mode approach to the perturba- 
tion equations will give eigenfunctions and eigenvalues 
but will not yield any other information, such as the sense 
of the induced flow, the nature of the transient motion, or 
the strength of the steady convecting state. 

Finite-amplitude nonlinear steady convection at supercrit- 
ical Rayleigh numbers has been addressed in the literature 
for many decades. A common approach is to use a Fourier 
series expansion to compute the final steady convective 
flow. The theoretical and experimental literature to 1980 
has been reviewed by Busse (1981). In particular, a few 
investigators attempted to track the evolution of discrete 
modes as the solution to an initial-value problem. 

Saltzman (1962) was apparently the first to solve the non- 
linear amplitude equations in an infinite horizontal layer 
and also to investigate the nature of the steady-state non- 
linear modes. Ogura and Yagihashi (1969), using a similar 
approach, found that initial conditions played a key role in 
determining the ultimate steady-state pattern. In a later 
paper Ogura (1971) investigated the effect of a variety of 
initial conditions (in the context of Fourier modes) on the 
final steady state. This study showed that the flow may 
evolve through a series of patterns before the final state is 
achieved and that the then-conventional rule that the final 
steady-state flow corresponded to the flow with maximum 
Nusselt number was false. 

In a recent and comprehensive work of this type, 
Goldhirsch et al. (1989) investigate natural convection in 
a rectangular enclosure by using Chebyshev expansions in 
the spatial coordinates and a time-split approach to the 
temporal coordinate. The flow is induced by a time- 
dependent exponential heating from below. The flow 
evolves by creating rolls at the boundaries that ultimately 
migrate into the cavity. 

In the present paper, a new finite-difference approach is 
used to solve the Boussinesq equations in primitive vari- 
ables. The method is based on a “discrete dispersion rela- 
tion” approach to the convective terms and a new splitting 
algorithm to construct compact implicit finite-difference 
formulas. The technique is used to examine the pattern 
selection process where finite-amplitude steady solutions 
are known to exist. Realistic boundary and initial condi- 
tions, consisting of insulated side walls and prescribed 
temperatures on the top and bottom walls, are imposed. A 
unit step in temperature along with “noisy” lower-wall 
boundary conditions are used to initiate the motion. The 
general flow features outlined by Goldhirsh et al. (1989) 
are confirmed by using the current finite-difference 
approach, but certain details of the instability growth dif- 
fer because of the use of insulated side walls, which 
induce a different set of eigenfunctions. 



The same finite-difference algorithm is used to correlate 
stability boundaries predicted by the classical perturbation 
method. The most (linearly) unstable mode is computed 
by treating the perturbed Boussinesq equations as a tran- 
sient stability problem. The eigenvalues and eigenfunc- 
tions, which appear as a natural outcome of the transient 
solution, are compared with direct solutions to the nonlin- 
ear equations. Numerical experiments in which the linear 
perturbation equations were used uncovered two instabil- 
ity modes over the range of Rayleigh numbers considered. 
Eigenfunctions from the linearized equations are also 
correlated with corresponding patterns computed during 
the nonlinear exponential growth and the finite-amplitude 
phases. The linearized solution is shown to dominate the 
flow evolution during the exponential growth phase, but it 
naturally fails to predict the ultimate nonlinear saturation. 

Numerically generated steady-state solutions confirm 
results from finite-amplitude theory and a parabolic bifur- 
cation emerges from the simulation. A series of computa- 
tions with noisy boundary conditions illustrates the role of 
transients in the pattern selection process. It is found that 
the final pattern is imprinted at an early stage of the flow 
where linearly unstable modes first become dominant. If 
two modes are possible, there is a chance that either may 
be activated. When a symmetric mode is generated it may 
be transformed to a stable rotational mode if both the 
amplitude and the symmetry of imposed disturbances are 
consistent with this pattern. Above the second critical 
Rayleigh number the three possible modes are clockwise 
or counterclockwise rotational modes and one symmetric 
mode consisting of two counterrotating vortices with a 
central updraft. 

Solutions to the Boussinesq equations are presented 1) as 
an example of a relatively complex fluid model to validate 
partially the algorithmic approach, 2) as a vehicle to 
examine some important physical effects, and 3) as a 
means of assessing the impact of linear theory on non- 
linear processes. Primitive variables are chosen so that 
fully three-dimensional problems can be computed with 
the same algorithm. In this report the equations of free 
convection are reviewed; the new finite-difference 
approach is described; numerical solutions are presented 
in terms of integrated quantities and instantaneous 
velocity-temperature fields; and finally, the process by 
which small disturbances are internalized is examined by 
using instantaneous snapshots of the flow during a critical 
time interval and comparing with a simple second-order 
dynamical systems model. The numerical algorithm used 
to simulate the flow is described in the appendix. 

The author would like to acknowledge the helpful discus- 
sions with Murray Tobak and Chris Hill regarding inter- 
pretation of the simulated flows. 
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convection speed 

Courant number, ax/h 

specific heat at constant pressure 

diffusion number 

length scale 

energy 

gravity vector 

space increment 

wavenumber in x-direction 

wavenumber in y-direction 

Prandtl number, v/% 

heat flux vector 

Rayleigh number, (JgATd 3 /v% 

critical Rayleigh number 

temperature 

time 

energy 

velocity vector 
eigenvalue (growth rate) 
coefficient of thermal expansion 
temperature increment 
kinematic viscosity 
density 

time increment 
radian frequency 
thermal diffusivity, k/p 0 c p 


Formulation of Equations and Solution 
Method 


Equations of Free Convective Flows 

The Boussinesq equations result from the following 
assumptions regarding the Navier Stokes and Energy 
equations for an incompressible fluid: 1) Fourier’s Law 
for an isotropic, homogeneous medium relates the heat 
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flux vector to the temperature gradient by q = -kVT; 

2) imposed temperature differences do not invoke large 
density changes; and 3) the thermodynamic state variables 
(density and energy) depend only on the temperature. 
Since the pressure is not considered a thermodynamic 
state variable, the pressure term in the equations of motion 
is considered a “hydrodynamic” quantity only. Density 
and energy changes induced by temperature fluctuations 
are computed from the equations p - po = -Ppo(T - To) 
and e - eo = c p (T - To), where the constants are the coef- 
ficient of thermal expansion (P) and the specific heat at 
constant pressure (c p ). 

With these assumptions, the equations of motion become 
the simpler set 

V • v* = 0 


Vv ' = 0 

— =-Vp , + PrV 2 7’ + RaPre , e (3) 

i^ = V 2 0' + 7 , e 

a 

The primed quantities in equation (3) are equivalent to 
those unprimed in equation (2) except for the temperature, 
which is related by 0 = 0' + (1 - y). These equations must 
be augmented with appropriate boundary conditions. The 
linear system represented by equation (3) is an eigenvalue 
problem and, unlike the full equations, it is solved with 
homogeneous boundary conditions. 


+ v* • Vv = -ZL_ + vV 2 v* + g - P(T* - To)g (1) 
at* po 

— +7* • VT* = -JS-VV 

at* pocp 

where starred symbols denote physical quantities. These 
equations are conventionally nondimensionalized with 
length scale d and time scale d 2 /x, where x = k/poc p is the 
thermal diffusivity. The final equations of motion in 
nondimensional form with the hydrostatic part eliminated 
are 


Numerical Method 

In two dimensions, equation (2) is a system of four partial 
differential equations for the dependent variables 
(u,v,0,p). The pressure plays a special role in the incom- 
pressible Navier-Stokes equations as the enforcer of mass 
conservation. If only those dependent variables with 
explicit time derivatives are retained, the following sys- 
tem of homogeneous matrix partial differential equations 
is obtained: 


au A au D au jfiv - n 

— + A — +B — +C +D — -+EU = 0 

0 t obc dy d \ 2 dy 2 


V • v = 0 

— + v • Vv = -Vp + PrV 2 v+RaPr0e (2) 

at 

— +v- ve=v 2 e 
at 

The velocity is scaled by x/d, the pressure by poWd) 2 . 
and the temperature change 0 by the driving temperature 
difference AT. The unit vector? is positive upward. Two 
key parameters appearing as coefficients are the Prandtl 
number (Pr = v/x) and the Rayleigh number 
(Ra as PgATd 3 /vx). In all the present computations it is 
assumed that the medium is air with Pr = 0.71. 


where U is a three-component solution vector consisting 
of (u,v,0), A = diag (u,u,u), B = diag(v,v,v), C = D = diag 
(_p r _Pr _1), and E = -Ra Pr 823. The first-order terms 
represent convection, the second-order terms diffusion, 
and the last is a body-force coupling the momentum and 
energy equations. The finite-difference algorithm is 
defined by examining the way that plane waves are pro- 
cessed by equation (4). The “ansatz” is a plane wave with 
arbitrary wave number k,l of the formU = Y(t) e ikx+il y, 
which is similar to that used in asymptotic wave theory 
(Whitham (1974)). Substituting U into equation (4) yields 
a first-order matrix ordinary differential equation for the 
amplitude function Y(t). 


The linear perturbation equations are obtained from equa- 
tion (2) by considering deviations from a basic state 
consisting of v = 0 and 0 = 1 - y, where y = 0 is the base 
of the cavity and 0 < (x,y) < 1. 


+ ikAY (t) + ilBY(t) - k 2 CY(t) - 1 2 DY(t) 

dt 

+ EY(t) = 0 (5) 

with a formal solution in terms of matrix exponentials: 
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Y(t) = YO exp(-J kAdt + ilBdt - k2Cdt- l 2 Ddt + Edt) 

( 6 ) 

The basis of the algorithm is to consider the change in Y 
over a small time interval x where temporal changes in the 
coefficient matrices are ignored. If a mesh parameter h 
and the ratio p = x/h is introduced, the integral is replaced 
by a product of matrix exponentials in terms of nondi- 
mensional spatial measures kh and lh. The wave parame- 
ter kh (or Ih) will play a fundamental role in the 
construction of finite -difference molecules. A significant 
challenge in deriving finite-difference algorithms is to 
maintain solution accuracy as kh (and Ih) varies over the 
allowable range (0 ,ji). Complete resolution is not possible 
as there will always be numerical artifacts leading to 
phase or amplitude distortion at short wavelengths. The 
goal of this effort is to derive formulas of maximum accu- 
racy for a given stencil in space-time. 


For convenience, the vectors Yo and Y are replaced by 
the conventional L? 1 and U 11 ** to represent the solution 
vector at sequential time steps. The solution is put in the 
special form 


u" +1 = e D^h 2 l 2 h 2 e Ci/h 2 k 2 h 2 e-EV iBp "V'Ap khu" (7) 


Several approximations have been made in the process of 
deriving this fundamental formula. First, the coefficient 
matrices A, B, etc. are assumed locally constant. Second, 
splitting the exponential into a product of terms is strictly 
valid for scalar arguments or matrices with special proper- 
ties. Both approximations would indicate a drop to first- 
order accuracy in the worst case. However, it is argued 
that an essential criterion of these simulations is to main- 
tain adequate phase resolution in the individual one- 
dimensional algorithms serving as the components of 
equation (7). 

Each term in equation (7) Is interpreted as an incremental 
one-dimensional operator pn U that gradually evolves 
the solution vector to U . The analogy with locally 
one-dimensional waves is complete when one considers 
the simple scalar wave operator du/cft + adu/dx with an 
input harmonic wave e ioM+ikx , In this case u n+1 = er* cn ** 
u n (from the dispersion relation © + ka = 0 and the 
Courant number cn = ax/h). Only one scalar operator is 
required to find u n+1 from u 11 . With reference to equa- 
tion (7), the matrix exponentials are replaced by partial 
solutions in a sequence of implicit finite-difference matrix 
operators: 

u 1 * 1 = (8) 


The solution sequence is a serial operation involving 
= C< x >U n , U® = etc. Each term has an 

obvious counterpart in equation (7) except for the inclu- 
sion of the pressure operator T . : The six-step process 
involves two convection sweeps, two diffusion sweeps, a 
velocity-temperature coupling, and a pressure 
computation. 

Numerical dispersion can introduce errors of a subtle 
nature; these errors are most apparent in wave problems. 
The nature of these effects and the associated waveform 
distortion effect using an acoustics model is discussed in 
Davis (1991). Numerical dispersion of nonlinear weakly 
dissipative systems is very difficult to isolate when all 
other interrelated physical effects are considered. How- 
ever, accurate modeling of the convective terms is critical 
to long-time solution accuracy. The procedure adopted is 
discussed in the previously cited reference where a dis- 
crete dispersion relation approach is used to obtain the 
best possible approximation for the chosen stencil. The 
ability of the algorithm to track accurately a range of 
wavelengths is also demonstrated in the cited reference. 
Boundary values for convection are computed from 
incoming and outgoing radiation conditions. It should be 
noted that these intermediate velocities do not satisfy the 
no-slip conditions. The complete convection algorithm is 
presented in the appendix. 

The coupling term is treated in a simple manner. Unlike 

the previous operators, only the local solution at (x,y) is 

required to advance to the next time level. Since no 

-*(3) 

derivatives are involved, the updated solution U is 
obtained from a Taylor series expansion of the matrix 
exponential e^ T to second order. The matrix E contains 
only one off-diagonal element, as indicated above. 

The pressure is introduced as a constraint to enforce mass 
conservation. Chorin (1968) first proposed an approach 
(which is now well accepted) to define a pressure term 
from the local velocity gradient by v^) - v@) = -x Vp so 
that p satisfies a Poisson equation V z p = V • 7x from 
the continuity condition V • v^ = 0. This divergence is 
computed using second-order central differences and the 
pressure computed from an available Poisson solver. 
Boundary conditions are imposed from the condition that 
the normal velocity v^) vanishes on all four walls. 

The diffusion operator is treated with an algorithm of the 
Crank-Nicolson form. Actually, the most accurate compu- 
tational molecule on a 3-point stencil is of order (t 2 , h^) in 
contrast to the Crank-Nicolson molecule, which is of 
order (t 2 , h 2 ). This algorithm also appears in the literature. 
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and it is referred to as the Douglas formula. (See Mitchell 
(1969).) However, the Crank-Nicolson formula was found 
to be adequate for the range of physical phenomena in this 
study. Diffusion is the last operation in the six-step pro- 
cess; it allows the imposition of no-slip boundary condi- 
tions. These algorithms are also shown in the appendix. 

The solution sequence is chosen in order to address the 
problem of intermediate boundary conditions in a physi- 
cally meaningful manner. The boundary conditions for 
convection are set by the operator itself. The coupling 
term causes no special boundary-condition problem since 
no neighboring mesh points are involved. The pressure 
plays a dual role of enforcing the conservation of mass 
and causing the normal component of boundary velocity 
to vanish. 

At this stage, an inviscid solution (the Euler equation) is 
available. The tangential velocity at the wall is a measure 
of the incremental wall vorticity strength at the indicated 
time step. The final diffusion sweeps define the viscosity- 
induced shear and enable all boundary conditions to be 
imposed exactly. During this final step the initially 
solenoidal velocity field remains solenoidal; this charac- 
teristic is a fundamental property of the heat equation. 

Solutions were obtained on a uniform 25 x 25 mesh with 
a time step chosen from x = 0.06h. There is no externally 
imposed velocity scale, but a Courant number Ui/h based 
on the final maximum convective velocity ranges to about 
one. These time steps are higher than usual since the con- 
vection algorithm is dissipation free and dispersion is low 
enough to minimize phase errors. The basic convection 
finite-difference molecule is phase accurate to a Courant 
number of about two, but dependent- variable coupling in 
the Boussinesq equations reduces its practical limit to a 
Courant number of about one. A typical 400-step run to 
steady flow takes about 90 seconds(s) on a Cray Y-MP. 


Numerical Solutions 

The initial-boundary-value problem is shown in figure 1 . 
The side walls are rigid and insulated to the flow of heat. 
The top and bottom walls are also impervious, but they 
are maintained at prescribed temperatures. Motion is 
initiated from rest by a unit step in temperature that causes 
a large energy input which induces a transient flow that 
evolves into a final steady-state pattern. This pattern can 
be either convecting or conducting, depending on 
Rayleigh number. This indicial approach is both 
physically plausible and firmly rooted in classical linear 
system theory. The temperature boundary condition on the 
lower wall is modeled with a “noisy heater;” that is, the 
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Figure 1. Geometry and typical flow pattern for natural 
convection in a box. 

temperature 1 + e(x,t) is considered to vary randomly in 
time (t > 0) with a prescribed space dependence. 

The response of a stable physical system governed by lin- 
ear homogeneous equations is an exponential decay to an 
equilibrium rest condition. The most significant parameter 
is its time constant. The situation for nonlinear systems is 
not as simple. The possibility of multiple steady solutions 
and bifurcation phenomena can easily be appreciated with 
the simple Ginzburg-Landau model considered by 
Rosenblat and Davis (1979). Solutions to the nonlinear 
model evolution equation 

jLLL = (R-4)U-U 3 

dt (9) 

U(0) = U 0 

are shown in figure 2(a) for a sequence of initial condi- 
tions Uo. If R = 2 (left panel), a value less than the critical 
R c = 4, the step response decays exponentially for any 
choice of initial conditions. When R = 6 (right panel), any 
disturbance ultimately resolves itself into the preferred 
steady solution U(«>) = ±1.414. If negative values of U are 
permitted, this is an example of a symmetric (pitchfork) 
bifurcation, the positive branch of which is depicted in 
figure 2(b). The branch U = 0 is unstable to infinitesimal 
perturbations if R > 4. 

The Rayleigh-Benard system can be interpreted in the 
spirit of the Ginzburg-Landau model with a global kinetic 
energy measure U(t) = jj 1 v vdA.IfRa<Rac, U(t) 
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(a) 





Figure 2. Response of a model nonlinear system to step inputs of 0.25, 0.50, 0.75, and 1.00. (a) System time history for 
subcritical (left) and supercritical (right) values of parameter R. (b) Bifurcation diagram of model nonlinear system. 


should decay exponentially whereas bifurcations will 
occur at supercritical Rayleigh numbers. The nature of 
Rayleigh-Benard instability is shown in figure 3, which 
depicts computations of U(t) on a semilog scale from 
numerical solutions of equation (2) with Ra as parameter. 
At the initial instant of motion the temperature step 
induces a uniform buoyant upflow that decays because of 
viscous stresses in the cavity and at the side walls. At a 
critical Rayleigh number buoyancy forces overwhelm 
dissipation causing exponential growth in U. At the 
critical Rayleigh number the growth is very slow and of 
long duration. There are four significant parameters 
relating to these curves: the onset time for growth, the 
growth rate, the time to nonlinear saturation, and the 
amplitude of the final state. 

This response is contrasted to linear perturbation theory 
that predicts the growth rate (eigenvalues) and some idea 
of the final state (the eigenfunctions). The linear solution 

is of the form e at f(x,y) where a is the eigenvalue and 
f(x,y) is its eigenfunction. This type of instability is 
known as absolute or global instability; it is illustrated in 


figure 4, which shows the time history of U from solutions 
of equation (3) as dashed lines. Three Rayleigh numbers 
are chosen and are superimposed on the related nonlinear 
solution from equation (2). Starting with a short transient, 
the curves representing solutions of the linear equations 
quickly assume exponential form and become unbounded 
in time. The exponential-growth portions of the curves are 
well predicted, but other significant parameters are left 
unresolved when linear theory is used. 1 

The growth rate a was extracted from the perturbation 
solutions, and it is shown in figure 5 by the curve labeled 
“mode 1 ” The interpolated value of Rayleigh number Ra 


J In the analysis of transition on airfoils and wings using linear 
stability theory an attempt to p redic t transition is often made 
based on the so-called “e-to-the-N-method ” This empirical 
procedure predicts the distance from instability to transition by 
assuming an amplitude growth ratio of e^ where N is about 10. 
This procedure is clearly analogous to the use of linear theory in 
figure 4 to predict the nonlinear saturation value. 
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Figure 3. Global behavior of Rayleigh-Benard system for 
Rayleigh numbers 2000, 3000, 4000, 5000, 7500, 
and 10,000. 



Figure 5. Curves of growth rate as a function of Rayleigh 
number for rotating (mode 1) and co-rotational (mode 2) 
eigenfunctions. Solid lines: 25 x 25 mesh, symbols: 
50x50 mesh. 



Time, I 

Figure 4 . Superposition of linear (- -) and nonlinear 

C ) response curves of the Rayleigh-Benard system 

for Rayleigh numbers 3000, 5000, and 10,000. 


at instability onset is 2512, which compares favorably 
with the values 2553 and 2585 reported by Lee et al. 
(1989) and Kurtzweg (1965), respectively. Numerical 
experiments with the perturbation equations revealed 
another instability mode that was only excited with per- 
turbations that were symmetric with respect to the line 
x = 0.5. The growth rates of this mode (mode 2 in the 
figure) are smaller than those of mode 1 and they become 
positive beyond a second critical Rayleigh number 
Ra^ =s 7017. Selected computations with a doubled grid 
are shown with filled symbols in figure 5. The eigenfunc- 
tions associated with each eigenvalue are obtained from 
the evolving numerical solution. During exponential 
growth they are simple multiples of one another. Figure 6 
shows computed linearized eigenfunctions for temper- 
ature and velocity corresponding to modes 1 and 2. (The 
temperature eigenfunctions are defined as deviations from 
a purely conducting linear state.) 

Steady-state solutions of the nonlinear equations may be 
represented on a bifurcation diagram such as that shown 
in figu re 7. Here, the square root of steady-state kinetic 
energy is taken as a representative amplitude. These 
equilibrium modes are computed by using initial condi- 
tions with the same symmetries as the linear eigen- 
functions. Modes 1 and 2 are al most exac t fits to the 
interpolating parabolas 0.1076 VRa-2512 and 
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Figure 7 . Finite-amplitude states for the Rayleigh-Benard 
system. 25x25 mesh, Ra = 7500. Symbols : computation, 
solid lines: best-fit parabolas. 


0.09176 YRa-7017. These results are consistent with 
Landau's analysis of nonlinear absolute stability (Landau 
and Lifshitz (1959), p. 104): they predicted that the ampli- 
tudes of the equilibrium state will be proportional to the 
square root of the deviation of Ra from the critical value. 
These results show that the numerical algorithm predicts 
both linear and nonlinear effects that are entirely consis- 
tent with the accepted behavior of such flows. 

In the next series of calculations the case Ra = 16,000 will 
be considered in some detail in order to distinguish the 
possible stable states associated with the linear eigenfunc- 
tions and the role of the boundary conditions in selecting 
each state. Two bifurcation curves were identified in fig- 
ure 7 as being representative of stable equilibrium states. 
The first is a unicellular vortex motion bifurcating from 
Ra = 2512 with an amplitude measure of 12.4, and the 
second a symmetric mode bifurcating from Ra = 7017 
with an amplitude measure of 8.7. 

A series of runs at Ra = 16,000 starting from rest with a 
small random 1 % temperature perturbation on the base of 
the cavity showed evidence of both modes. The vortex 
modes (either clockwise or counterclockwise rotation) 
appeared more than twice as often as the two-vortex con- 
figuration. Energy growth curves (on a linear scale) are 
shown in figure 8. There seemed to be some indecision 
just after the exponential-growth phase as to which mode 
would predominate. In fact, this case shows that the 



Time, t 

Figure 8. Time evolution of global parameter U showing 
selection of mode 1 and mode 2 final states. 

25x25 mesh, Ra = 16,000. 

mode 1 solution coincides with the mode 2 growth-rate 
path at early times. The final state for both modes is pre- 
sented in figure 9, with panels of perturbation tempera- 
ture, total temperature, and velocity vectors. These modes 
retain remnants of the linear eigenfunctions of figure 6, 
but they are distorted by the nonlinear convection term in 
the Navier-Stokes equation. 

Examination of the transient motion shows that the non- 
linear stability curves possess a characteristic energy dip 
(cf. fig. 3) where many significant events occur. In this 
region the flow first aligns itself into the preferred spatial 
pattern (an “immature eigenfunction”) in preparation for 
exponential growth. The field of velocity vectors is exam- 
ined at the six points indicated at the top of figure 10. The 
first three represent the internalization period of the initial 
conditions (also called “receptivity”), the fourth the 
minimal energy state, and the last two the instability mode 
that determines the final state. 

A generally uniform updraft in figure 10(a) is the initial 
response to the temperature step function; the large shear 
generates side-wall vorticity. The velocity terms in the 
Navier-Stokes equation convect and diffuse this vorticity, 
but the pressure plays a special role. In enforcing continu- 
ity, it tends to turn or rotate the flow. The net effect, along 
with the no-slip boundary conditions, is to generate two 
narrow recirculating eddies as shown in figure 10(b). The 
vortices migrate from the wall and align themselves more 
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Figure 9. Maps of mode 1 and 2 final states. 25 x 25 mesh, Ra = 16,000. 
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or less symmetrically in the manner of a mode 2 eigen- 
function as in figure 10(c). This spatial realignment takes 
place in a field of decreasing energy. (All velocities in 
fig. 10 are to the same scale.) The flow is now configured 
for unstable growth. The small, random, 1% temperature 
perturbations on the lower boundary buffet the eddies 
about, and a chance push in the right direction induces 
exponential growth. Panel 10(d) shows the point at which 
these disturbances first start to affect the flow and pan- 
els 10(e) and 10(0 confirm the dominance of the left vor- 
tex. This counterclockwise eddy soon extinguishes its 
partner and the final state will consist of a single large 
vortex. 

Some investigations into what we have called the “sym- 
metric state” are shown in figures 1 1 and 12. Unlike 
previous examples, the flow is first established as the 
symmetric two-vortex mode before disturbances are 
imposed. Disturbances are initiated at t = 0.5 and released 
at t = 0.6. Three types of temperature distributions along 
the bottom wall are chosen: 

(a) 0(x) = eRand{-l,l) (x -0.5) 

(b) 0(x) = £ (10) 

(c) 0(x) = e (x - 0.5) 

where Rand{-1 , 1 } is a random number in the indicated 
range. The disturbances (a) have no spatial coherence 
while disturbances (b) and (c) are symmetric and antisym- 
metric, respectively, with respect to the vertical line of 
symmetry x = 0.5. Figure 1 1 indicates that random pertur- 
bations with peak-to-peak amplitudes as much as 100% 
cause the two vortices to react quite violently with one 
another but do not change the fundamental modal struc- 
ture, which settles down to its initial configuration. 
Symmetric disturbances of type (b) cause violent excur- 
sions, but they also fail to upset the basic symmetry. 

Case (c) of equation (10), presented in figure 12, shows 
that above a critical amplitude the flow becomes asym- 
metric and evolves into the stable vortex mode state. In 
summary, 1) the symmetric state is stable to symmetric 
disturbances of any size, 2) the symmetric state is stable 
to small antisymmetric or random disturbances, and 3) the 
symmetric state is unstable to larger antisymmetric 
disturbances. 

The mode selection process will now be investigated with 
the aid of a model second-order nonlinear dynamic sys- 
tem. Consider a phase-plane (or state-space) representa- 
tion with the kinetic energy U(t) as ordinate and dU/dt as 
abscissa. The instantaneous state actually requires an 
infinite-dimensional velocity vector, but U will represent 
a convenient measure for a model second-order damped 
nonlinear oscillator. Finite-amplitude steady flows consist 

of discrete limit points on the U = 0 axis. At supercritical 



Figure 1 1 . Energy time history for random temperature 
inputs to base of cavity during time interval 0.5 - 0.6 with 
various amplitudes. 



Time, t 

Figure 12. Energy time history for anti-symmetric 
temperature inputs to base of cavity during the time 
interval 0.5 - 0.6 with various amplitudes. Ra = 16,000. 
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Rayleigh numbers there may be many limit points: singu- 
larities possessing domains of attraction to which all solu- 
tion paths congregate. Figure 12 shows that there are 
attractors at U = 76 and U = 154 for Ra = 16,000. The 
state-space path corresponding to £ = 1 is shown in 
Figure 13 to be attracted from rest to the mode 2 state at A 
and then perturbed to a path that encroaches on the 
mode 1 domain of attraction. Points A and B along the 
path correspond to the imposition and elimination of the 
disturbance at t = 0.5 and t = 0.6. The curve beyond 
point B is the response of the Boussinesq equations to 
initial conditions at B that were vigorous enough to propel 
the solution path away from the mode 2 attractor. 

The stable state that will be realized depends on the value 
of the parameter e; there must be a critical amplitude 
depending on £ that separates these flows. The particular 
path corresponding to this critical value should pass 
through a saddle-type singular point in the phase plane. 
Numerical experiments using the bisection method identi- 
fied critical amplitude bounds in the range 
0.9031 < 8 < 0.9047. State-space trajectories for four 



Rate of energy, U(t) 


Figure 13. Phase-plane representation for flow evolving 
from rest to a mode 2 symmetric state at U= 76 and 
perturbed to a mode 1 state at U = 154 with anti-symmetric 
disturbances. Equation 10(c), e= 1.0. The path moves 
below an unstable singular point at C. 


amplitudes in the vicinity of the critical value are shown 
in figure 14. The trajectories pass near a singularity at 
point C, representing an unstable equilibrium. The role of 
this unstable state is appreciated from the time histories in 
figure 15 for £ = 0.90312 and £ = 0.9047. The system 
responds along a common path after removal of the 
perturbation at point B until it encounters the instability at 
U = 53.6; from here they evolve to a final state fated by 
the disturbance history. A change of 0.17% in £ induces a 
final state with an energy change of over 200%. A flow 
pattern near the unstable equilibrium point is shown in 
figure 16. Here the velocities prefer a diagonally sym- 
metric pattern with vortices of opposite sense. The pattern 
for a sign-reversed flow near the instability point is 
similar except that vortices appear at opposite comers 
with sign-reversed flow. 

The singular points of the trajectories shown in figures 13 
and 14 do not line up in the usual stable-unstable-stable 
sequence because of the use of the positive-definite quan- 
tity U as the characteristic measure. A more conventional 
diagram can be recovered if a dependent variable that can 
be either positive or negative is selected as the character- 
istic measure. Figure 17 shows transient solution paths to 
the static attractors associated with the quantity A0 (the 
temperature difference across the box at midheight). The 
sense of the single vortex solution can be positive or 
negative, as revealed by the symmetric location of the 
spiral attractors. The origin of this phase diagram repre- 
sents both the initial conducting flow and the symmetric 
mode 2 state. The boundaries of the basins of attraction 
are defined by separatrices; they are approximately indi- 
cated by the solid lines, which represent trajectories for 
£ = ±0.9047. Representative paths in the two basins of 
attraction for positive A0 are shown by dashed lines. 

This digression to nonlinear dynamics for a specific dis- 
turbance mode and Rayleigh number is presented to illus- 
trate how the numerical approach can be used to examine 
and interpret parametric trends with a global measure. A 
complete analysis of all possible modes, phase -plane 
parameters, and flow states is beyond the scope of this 
paper. 

The examples discussed above for Ra > Rac C have a 
richer structure than the response for Ra < Ra^ that 
allows only simple rotational-flow solutions. If the 
Rayleigh number is greater than about 20,000, fine-scale 
secondary instabilities appear and the role of three- 
dimensional disturbances and nonlinear oscillatory modes 
must be considered ab initio. 
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Figure 14. Phase-plane solution paths for disturbances that pass close to an instability point at C. (a) e = 0.9000, 
(b) e = 0.90312, (c) £ = 0.90469, (d) £ = 0.91250. 
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Figure 15. Energy time histones for amplitudes that pass 
near the unstable equilibrium. 


Velocity vectors 



Figure 16. Velocity-vector maps in vicinity of unstable 
equilibrium att= 0.6. 



A ®mid 

Figure 17. Phase-plane trajectories for solution paths 
using temperature differences across box at midheight. 


Conclusions 

A new primitive-variable finite-difference simulation of 
the Boussinesq equations has revealed details of the 
initial-flow structure and has traced modes of instability 
from inception to the final nonlinear state. The finite- 
difference scheme used in these computations differs from 
other techniques in the use of matrix-exponential splitting 
to separate physical and dimensional processes. This 
method allows construction of efficient, low-dispersion 
computational molecules as well as relatively coarse 
meshes and longer time steps. 

As a result of this study, a picture emerges of the response 
of the Boussinesq model to imposed disturbances. The ini- 
tial response involves a realignment of the velocity field 
that is strongly influenced by the incompressibility con- 
straint. Once the field attains a small but spatially coher- 
ent pattern, strong coupling with other dependent 
variables induces a rapid and large-scale change in the 
flow energy. Nonlinear effects cause saturation to a steady 
flow governed by the magnitude of the externally applied 
temperature field. While the linear perturbation equations 
can accurately predict the growth rates and eigenfunction 
patterns, the nonlinear equations must be analyzed in 
detail before the initial internalization process and the 
process of final pattern selection can be understood. 
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In this study only step-type motions — specifically instan- 
taneously applied heating — have been used. Other initial 
conditions may have different initial responses, but the 
nature of instability growth will not change. In addition, 
only two-dimensional flows were considered. Expanding 
the parameter space to a third dimension and to other 
aspect ratios may invoke other interesting modes. The 
Boussinesq model of flow stability couples the dependent 
variables in a particular but relatively simple manner. 
Other models, especially those governing convective 
flows with inflow-outflow boundaries, will have other 
modal responses which will be investigated in the future 
using the numerical techniques introduced. 

Appendix 

Natural Convection Algorithm 

The algorithm used to simulate the Boussinesq equations 
is presented schematically in equation (7). The procedure 
is similar to other time-splitting approaches except that 
each operator and each coordinate direction is identified 
as a separate physical process. It has been difficult to 
obtain accurate composite finite-difference formulas for 
convection and diffusion. However, simpler algorithms 
are obtained by separating the formulas into distinct steps. 
In this appendix the computational molecule for each step 
is shown in some detail. 

The convection algorithm is split into separate x- and 
y-sweeps. The system is closed with finite-difference 
approximations to radiation conditions; this approach 
avoids the question of intermediate boundary conditions. 
Consider, for example, a uniform mesh at time t = nx and 
points x = (i - l)h, y = (j - l)h (1 < i, j < N) with velocity 
denoted by uj 1 -. A local Courant number is defined in 
terms of the local velocity by cn = u^-x/h. The algorithm 
for updating u at an interior mesh point during an x-sweep 
is: 

-J- (cn - l)(cn - 2)u?Vi ~ L (cn - 2) 

12 ‘‘ 1J 6 

x (cn + 2)ii Pt 1 +-L (cn + l)(cn + 2)^ 

= -L(cn + l)(cn + 2)i^ . ; (Al) 

12 1 d 

-l(cn-2)(cn + 2)u?j 

+ kF ( Cn - l^n - 2)*^+l j 

The boundary algorithms at i = 1 and i = N are 


(cn - l)Uj j 1 - (cn + Dujj 1 = 

-(cn + l)Uj • + (cn- l)u2j 

d d (A2) 

-(cn - l)ujjjt| j + (cn + = 

(cn+ 1 )u^f_ i j -(cn - l)u|Jjj 

These sweeps convect the dependent variable in x without 
specific reference to the boundary conditions. If an arbi- 
trary plane wave with wavenumber k is propagated 
through the mesh, equation (Al) is accurate to (kh) 4 . 
Equation (A2) involves only two mesh points and it is 
accurate to (kh)^. In the commonly used terminology, 
these equations are accurate in space-time to fourth- and 
second-order, respectively. This criterion is not a com- 
plete representation of the accuracy of the finite- 
difference approximation. In actual simulations the 
accuracy attainable with a given mesh depends on the 
gradients of the computed solution. In any physical 
problem a higher accuracy algorithm will always be a 
better choice, but it is more difficult to decide whether or 
not the computed solution is a more valid physical 
representation. The algorithm for v and 0 is exactly the 
same as the preceding equations. 

The y-sweep is computed in the same manner as equa- 
tions (Al) and (A2) except that cn is now defined by 
vf-x/h with the j index being variable. The notation 
(n + 1) does not mean that the solution bears any 
resemblance to the actual physical solution at time n + 1; 
it only denotes a partial correction as one step toward the 
ultimate physical solution. 

The vector of dependent variables is now corrected for the 
coupling term involving the matrix E in equation (7). This 
term, which does not involve any spatial derivatives, is 
equivalent to the simple matrix equation dltfdl + EU = 0. 
The updated vector is obtained from a Taylor series in 
time with X as the small parameter: 

U = Uoe-Ex = Uo(I + ET + lEEx 2 + ...) (A3) 

2 

The matrix E is very simple; all powers of E vanish 
except for the first. The vertical velocity is updated from 
the simple formula 

v !j' 1 = v i I j(l +E 2 3'0 (A4) 

Note that this revised velocity does not involve the solu- 
tion at neighboring mesh points. 

Updated velocities u, v are now available, and the pres- 
sure is computed in order that the velocities remain com- 
patible with the incompressibility condition. The Poisson 
equation for p as discussed in the body of the paper is: 
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V 2 p = J^9u/3x + dv/^y) 

— = u, (x = 0, x = 1, 0 <y < 1) 
dx 

|P = v,(y = 0,y=l,0^x<l) 
dy 


where the diffusion number dn = Prt/h^ for velocity dif- 
fusion and dn = i/h 2 for temperature diffusion. The 
(A5) formula is exactly the same for y-diffusion. Boundary 

conditions based on the no-slip condition for velocity and 
either the no-slip or the no-flux condition for temperature 
are used. 


This two-dimensional elliptic boundary-value problem on 
a uniform mesh is solved with available software library 
routines. The solution for p is used to update the velocities 
from the defining equation 


u new = -idp/dx + u old 
v new = — /dy + v 0 ld 


(A6) 


The last step is a diffusion to the final velocity as deter- 
mined from a Crank-Nicolson sweep. In the x-direction: 

-Idnu^+a.dntu^'-ldnu?^ 

=2 d "u"-lj +a - dn Kj*i- dnu "*lj <A7) 


Each of the tridiagonal equations appearing in these equa- 
tions is solved by using variants of the Thomas algorithm. 
Although the diffusion equation is diagonally dominant 
and causes no solution difficulties, the Thomas algorithm 
will fail for equation (Al) if some of the local values of cn 
exceed unity. This failure is not due to the algorithm, 
which was shown in the cited references to be accurate to 
at least cn = 2, but it is caused by a breakdown in diagonal 
dominance required by the algorithmic construction. This 
breakdown forces the allowable time step to be reduced 
somewhat. Other classes of sparse matrix solvers that 
might alleviate this restriction are under investigation. 
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